In this report, we reproduce the analyses in the follow-up behavioral study 2.

prep data

First, we load the relevant packages, define functions and plotting aesthetics, and load and tidy the data.

load packages

if(!require('pacman')) {
    install.packages('pacman')
}

pacman::p_load(tidyverse, purrr, fs, knitr, lmerTest, ggeffects, parameters, kableExtra, boot, devtools, EMAtools, install = TRUE)
devtools::install_github("hadley/emo")

define functions

# MLM results table function
table_model = function(model_data, eff_size = FALSE, word_count = TRUE) {
  
  results = model_data %>%
    broom.mixed::tidy(conf.int = TRUE) %>%
    filter(effect == "fixed") %>%
    rename("SE" = std.error,
           "t" = statistic,
           "p" = p.value) %>%
    select(-group, -effect) %>%
    mutate_at(vars(-contains("term"), -contains("p")), round, 2) %>%
    mutate(term = gsub("article_cond", "", term),
           term = gsub("\\(Intercept\\)", "control", term),
           term = gsub("sharing_type", "sharing type", term),
           term = gsub("msg_rel_self_z", "self-relevance", term),
           term = gsub("msg_rel_social_z", "social relevance", term),
           term = gsub("conditiontimed", "group (timed)", term),
           term = gsub("conditionuntimed", "group (untimed)", term),
           term = gsub("contentclimate", "content (climate)", term),
           term = gsub("siteUSA", "sample (USA)", term),
           term = gsub("n_c", "word count", term),
           term = gsub(":", " x ", term),
           p = ifelse(p < .001, "< .001",
               ifelse(p > .999, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))),
           `b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", estimate, conf.low, conf.high)) 
  
  if (word_count == TRUE) {
    results = results %>%
      mutate(term = gsub("control", "intercept", term))
  }
  
  if (eff_size == TRUE) {
    eff_size = lme.dscore(model_data, data = data, type = "lme4") %>%
      rownames_to_column(var = "term") %>%
      mutate(term = gsub("article_cond", "", term),
             term = gsub("article_cond", "", term),
             term = gsub("\\(Intercept\\)", "control", term),
             term = gsub("sharing_type", "sharing type", term),
             term = gsub("msg_rel_self_between", "self-relevance", term),
             term = gsub("msg_rel_social_between", "social relevance", term),
             term = gsub("contentclimate", "content (climate)", term),
             term = gsub(":", " x ", term),
             d = sprintf("%.2f", d)) %>%
      select(term, d)
    
    results %>%
      left_join(., eff_size) %>%
      mutate(d = ifelse(is.na(d), "--", d)) %>%
      select(term, `b [95% CI]`, d, df, t, p)
    
  } else {
    results %>%
      select(term, `b [95% CI]`, df, t, p)
  }
}

# simple effects function
simple_effects = function(model, sharing = FALSE) {
  if(sharing == FALSE) {
    results = emmeans::contrast(emmeans::emmeans(model, ~ article_cond | condition),
                            "revpairwise", by = "condition", adjust = "none") %>%
      data.frame() %>%
      filter(grepl("control", contrast)) %>%
      rename("group" = condition) %>%
      select(contrast, group, estimate, p.value)
  } else {
    results = emmeans::contrast(emmeans::emmeans(model, ~ article_cond | condition + sharing_type),
                            "revpairwise", by = "condition", adjust = "none") %>%
      data.frame() %>%
      filter(grepl("- control", contrast)) %>%
      filter(!grepl("^control", contrast)) %>%
      extract(contrast, c("exp_sharing", "control_sharing"), ".* (0|1) - control (0|1)", remove = FALSE) %>%
      filter(exp_sharing == control_sharing) %>%
      mutate(sharing_type = ifelse(exp_sharing == 0, "broadcast", "narrowcast"),
             contrast = gsub("0|1", "", contrast)) %>%
      rename("group" = condition) %>%
      select(contrast, sharing_type, group, estimate, p.value)
  }
  
  results %>%
    mutate(p.value = ifelse(p.value < .001, "< .001",
                      ifelse(p.value == 1, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p.value))))) %>%
    kable(digits = 2) %>%
    kableExtra::kable_styling()
}

define aesthetics

palette_condition = c("self" = "#ee9b00",
                      "control" = "#bb3e03",
                      "other" = "#005f73")
palette_dv = c("self-relevance" = "#ee9b00",
               "social relevance" = "#005f73",
               "sharing" = "#56282D")

plot_aes = theme_minimal() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        text = element_text(size = 16, family = "Futura Medium"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(color = "black"),
        axis.line = element_line(colour = "black"),
        axis.ticks.y = element_blank())

load data

data = read.csv("../data/study2_data.csv", stringsAsFactors = FALSE) %>%
  mutate(article_cond = ifelse(article_cond == "social", "other", article_cond),
         msg_rel_self_z = scale(msg_rel_self, center = TRUE, scale = TRUE),
         msg_rel_social_z = scale(msg_rel_social, center = TRUE, scale = TRUE)) %>%
  rename("condition" = group) %>%
  filter(sharing_type_key == "msg_sharing_narrow")

sub_conditions = data %>%
  select(SID, condition) %>%
  unique()

descriptives

group ns

Sample size by group

data %>%
  select(condition, SID) %>%
  unique() %>%
  group_by(condition) %>%
  summarize(n = n()) %>%
  kable() %>%
  kable_styling()
condition n
comment 124
timed 157
untimed 167

ratings

Summarize means and SDs

data %>%
  gather(variable, value, msg_share, msg_rel_self, msg_rel_social) %>%
  group_by(variable) %>%
  summarize(M = mean(value, na.rm = TRUE),
            SD = sd(value, na.rm = TRUE)) %>%
  mutate(variable = ifelse(variable == "msg_rel_self", "self-relevance",
                    ifelse(variable == "msg_rel_social", "social relevance", "sharing intention"))) %>%
  kable(digits = 2) %>%
  kableExtra::kable_styling()
variable M SD
self-relevance 49.33 34.19
social relevance 54.47 33.04
sharing intention 35.65 33.96

H2

Do the manipulations increase relevance? Is this effect stronger in the comment group?

self-relevance

✅ H2a: Self-focused intervention (compared to control) will increase self-relevance

We replicate our previous work in the comment group: the self-focused condition increases self-relevance compared to the control

✅ This effect is smaller for both the timed and untimed groups

mod_h2a = lmer(msg_rel_self ~ 1 + article_cond * condition + (1 + article_cond | SID),
              data = data,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h2a = table_model(mod_h2a)

table_h2a %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 49.91 [45.62, 54.20] 456.24 22.85 < .001
other 4.03 [0.43, 7.62] 462.82 2.20 .028
self 14.44 [10.47, 18.40] 457.92 7.15 < .001
group (timed) -3.52 [-9.23, 2.19] 449.41 -1.21 .227
group (untimed) -4.66 [-10.31, 0.98] 450.11 -1.62 .105
other x group (timed) -5.68 [-10.44, -0.92] 451.05 -2.35 .019
self x group (timed) -9.83 [-15.10, -4.56] 448.78 -3.67 < .001
other x group (untimed) -3.26 [-7.96, 1.45] 451.42 -1.36 .174
self x group (untimed) -11.66 [-16.86, -6.46] 449.43 -4.40 < .001

simple effects by group

simple_effects(mod_h2a)
contrast group estimate p.value
other - control comment 4.03 .028
self - control comment 14.44 < .001
other - control timed -1.66 .296
self - control timed 4.61 .009
other - control untimed 0.77 .617
self - control untimed 2.77 .105

summary

summary(mod_h2a)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_self ~ 1 + article_cond * condition + (1 + article_cond |  
##     SID)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 51167.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.90817 -0.68520  0.05154  0.68778  2.92000 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr       
##  SID      (Intercept)       392.42   19.810              
##           article_condother  25.35    5.035   -0.12      
##           article_condself  117.09   10.821   -0.20  0.52
##  Residual                   739.97   27.202              
## Number of obs: 5312, groups:  SID, 448
## 
## Fixed effects:
##                                    Estimate Std. Error      df t value
## (Intercept)                          49.907      2.184 456.241  22.852
## article_condother                     4.026      1.829 462.817   2.201
## article_condself                     14.436      2.019 457.922   7.149
## conditiontimed                       -3.519      2.906 449.406  -1.211
## conditionuntimed                     -4.663      2.871 450.111  -1.624
## article_condother:conditiontimed     -5.684      2.422 451.052  -2.347
## article_condself:conditiontimed      -9.829      2.680 448.778  -3.668
## article_condother:conditionuntimed   -3.255      2.393 451.421  -1.360
## article_condself:conditionuntimed   -11.661      2.647 449.434  -4.405
##                                                Pr(>|t|)    
## (Intercept)                        < 0.0000000000000002 ***
## article_condother                              0.028245 *  
## article_condself                       0.00000000000349 ***
## conditiontimed                                 0.226689    
## conditionuntimed                               0.105002    
## article_condother:conditiontimed               0.019358 *  
## article_condself:conditiontimed                0.000274 ***
## article_condother:conditionuntimed             0.174483    
## article_condself:conditionuntimed      0.00001325096630 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                     (Intr) artcl_cndt artcl_cnds cndtnt cndtnn
## artcl_cndth         -0.425                                    
## artcl_cndsl         -0.442  0.494                             
## conditintmd         -0.751  0.319      0.332                  
## conditnntmd         -0.761  0.323      0.336      0.572       
## artcl_cndthr:cndtnt  0.321 -0.755     -0.373     -0.420 -0.244
## artcl_cndslf:cndtnt  0.333 -0.372     -0.754     -0.437 -0.253
## artcl_cndthr:cndtnn  0.325 -0.764     -0.377     -0.244 -0.420
## artcl_cndslf:cndtnn  0.337 -0.377     -0.763     -0.253 -0.438
##                     artcl_cndthr:cndtnt artcl_cndslf:cndtnt artcl_cndthr:cndtnn
## artcl_cndth                                                                    
## artcl_cndsl                                                                    
## conditintmd                                                                    
## conditnntmd                                                                    
## artcl_cndthr:cndtnt                                                            
## artcl_cndslf:cndtnt  0.490                                                     
## artcl_cndthr:cndtnn  0.577               0.284                                 
## artcl_cndslf:cndtnn  0.285               0.575               0.491

social relevance

✅ H2b: Other-focused intervention (compared to control) will increase social relevance

We replicate our previous work in the comment group: the other-focused condition increases social relevance compared to the control

✅ This effect is smaller for both the timed and untimed groups

mod_h2b = lmer(msg_rel_social ~ 1 + article_cond * condition + (1 + article_cond | SID),
              data = data,
              control = lmerControl(optimizer = "bobyqa"))

model table

table_h2b = table_model(mod_h2b)

table_h2b %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 53.76 [49.40, 58.12] 455.12 24.22 < .001
other 14.45 [10.91, 17.99] 459.95 8.03 < .001
self 9.87 [6.28, 13.47] 459.93 5.39 < .001
group (timed) -3.11 [-8.92, 2.70] 449.05 -1.05 .293
group (untimed) -4.36 [-10.10, 1.38] 449.79 -1.49 .136
other x group (timed) -13.98 [-18.67, -9.29] 448.90 -5.86 < .001
self x group (timed) -6.86 [-11.63, -2.09] 450.70 -2.83 .005
other x group (untimed) -10.26 [-14.89, -5.62] 449.24 -4.35 < .001
self x group (untimed) -7.39 [-12.11, -2.68] 451.27 -3.08 .002

simple effects by group

simple_effects(mod_h2b)
contrast group estimate p.value
other - control comment 14.45 < .001
self - control comment 9.87 < .001
other - control timed 0.47 .765
self - control timed 3.01 .059
other - control untimed 4.19 .006
self - control untimed 2.48 .110

summary

summary(mod_h2b)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_social ~ 1 + article_cond * condition + (1 + article_cond |  
##     SID)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 50361.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1317 -0.5822  0.0851  0.6392  3.0888 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr       
##  SID      (Intercept)       444.83   21.091              
##           article_condother  76.97    8.773   -0.14      
##           article_condself   91.86    9.585   -0.26  0.47
##  Residual                   615.48   24.809              
## Number of obs: 5312, groups:  SID, 448
## 
## Fixed effects:
##                                    Estimate Std. Error      df t value
## (Intercept)                          53.758      2.220 455.116  24.217
## article_condother                    14.450      1.801 459.953   8.026
## article_condself                      9.875      1.830 459.932   5.395
## conditiontimed                       -3.111      2.956 449.046  -1.052
## conditionuntimed                     -4.362      2.920 449.786  -1.494
## article_condother:conditiontimed    -13.982      2.386 448.900  -5.860
## article_condself:conditiontimed      -6.861      2.428 450.703  -2.825
## article_condother:conditionuntimed  -10.256      2.358 449.245  -4.350
## article_condself:conditionuntimed    -7.394      2.399 451.274  -3.082
##                                                Pr(>|t|)    
## (Intercept)                        < 0.0000000000000002 ***
## article_condother                    0.0000000000000085 ***
## article_condself                     0.0000001100232017 ***
## conditiontimed                                  0.29324    
## conditionuntimed                                0.13588    
## article_condother:conditiontimed     0.0000000089507683 ***
## article_condself:conditiontimed                 0.00493 ** 
## article_condother:conditionuntimed   0.0000168795897815 ***
## article_condself:conditionuntimed               0.00218 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                     (Intr) artcl_cndt artcl_cnds cndtnt cndtnn
## artcl_cndth         -0.388                                    
## artcl_cndsl         -0.433  0.500                             
## conditintmd         -0.751  0.291      0.325                  
## conditnntmd         -0.760  0.295      0.329      0.571       
## artcl_cndthr:cndtnt  0.293 -0.755     -0.378     -0.383 -0.222
## artcl_cndslf:cndtnt  0.326 -0.377     -0.754     -0.428 -0.248
## artcl_cndthr:cndtnn  0.296 -0.764     -0.382     -0.222 -0.383
## artcl_cndslf:cndtnn  0.330 -0.382     -0.763     -0.248 -0.429
##                     artcl_cndthr:cndtnt artcl_cndslf:cndtnt artcl_cndthr:cndtnn
## artcl_cndth                                                                    
## artcl_cndsl                                                                    
## conditintmd                                                                    
## conditnntmd                                                                    
## artcl_cndthr:cndtnt                                                            
## artcl_cndslf:cndtnt  0.497                                                     
## artcl_cndthr:cndtnn  0.576               0.288                                 
## artcl_cndslf:cndtnn  0.288               0.575               0.497

combined plot

predicted_h2a = ggeffects::ggpredict(mod_h2a, c("article_cond", "condition")) %>%
              data.frame() %>%
  mutate(model = "self-relevance")

predicted_h2b = ggeffects::ggpredict(mod_h2b, c("article_cond", "condition")) %>%
              data.frame() %>%
  mutate(model = "social relevance")

bind_rows(predicted_h2a, predicted_h2b) %>%
  mutate(x = factor(x, levels = c("self", "control", "other"))) %>%
  ggplot(aes(x = group, y = predicted, color = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1) +
  facet_grid(~ model) +
  coord_flip() +
  scale_color_manual(name = "", values = palette_condition) +
  labs(x = "", y = "\nmean predicted relevance rating") +
  plot_aes +
  theme(legend.position = "top")

H3

Is greater self and social relevance associated with higher sharing intentions?

✅ H1a: Greater self-relevance ratings will be associated with higher narrowcast sharing intentions

✅ H1a: Greater social relevance ratings will be associated with higher narrowcast sharing intentions

mod_h3 = lmer(msg_share ~ msg_rel_self_z + msg_rel_social_z + (1 + msg_rel_self_z + msg_rel_social_z | SID),
              data = data,
              control = lmerControl(optimizer = "bobyqa"))

plot

vals = seq(-2,2,.1)

predicted_h3 = ggeffects::ggpredict(mod_h3, c("msg_rel_self_z [vals]")) %>%
  data.frame() %>%
  mutate(variable = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h3, c("msg_rel_social_z [vals]")) %>%
              data.frame() %>%
              mutate(variable = "social relevance"))

predicted_sub_h3 = ggeffects::ggpredict(mod_h3, terms = c("msg_rel_self_z [vals]", "SID"), type = "random") %>%
  data.frame() %>%
  mutate(variable = "self-relevance") %>%
  bind_rows(ggeffects::ggpredict(mod_h3, c("msg_rel_social_z [vals]", "SID"), type = "random") %>%
              data.frame() %>%
              mutate(variable = "social relevance"))

predicted_h3 %>%
  ggplot(aes(x, predicted)) +
  stat_smooth(data = predicted_sub_h3, aes(group = group, color = variable),
              geom ='line', method = "lm", alpha = .05, linewidth = .75, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = variable), alpha = .5, color = NA) +
  geom_line(aes(color = variable), size = 1.5) +
  facet_grid(~variable) +
  scale_color_manual(name = "", values = palette_dv) +
  scale_fill_manual(name = "", values = palette_dv) +
  labs(x = "\nrelevance rating", y = "predicted sharing intention rating\n") +
  plot_aes +
    theme(legend.position = "none")

model table

table_model(mod_h3) %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 33.88 [32.15, 35.60] 436.70 38.63 < .001
self-relevance 3.42 [2.23, 4.62] 294.46 5.63 < .001
social relevance 16.55 [15.23, 17.87] 358.58 24.71 < .001

summary

summary(mod_h3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_share ~ msg_rel_self_z + msg_rel_social_z + (1 + msg_rel_self_z +  
##     msg_rel_social_z | SID)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 46155.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0902 -0.3691 -0.0377  0.3429  5.5329 
## 
## Random effects:
##  Groups   Name             Variance Std.Dev. Corr       
##  SID      (Intercept)      309.09   17.581              
##           msg_rel_self_z    62.52    7.907    0.27      
##           msg_rel_social_z  94.80    9.736    0.76 -0.39
##  Residual                  262.07   16.189              
## Number of obs: 5312, groups:  SID, 448
## 
## Fixed effects:
##                  Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)       33.8769     0.8770 436.6954  38.629 < 0.0000000000000002 ***
## msg_rel_self_z     3.4207     0.6072 294.4623   5.634         0.0000000412 ***
## msg_rel_social_z  16.5515     0.6699 358.5845  24.707 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) msg_rl_sl_
## msg_rl_slf_  0.161           
## msg_rl_scl_  0.536 -0.617

H5

Do the manipulations increase sharing intentions? Is this effect stronger in the comment group?

Here we focus on narrowcasting only since that is the type of sharing we used in fMRI study 1.

✅ H5a: Self-focused intervention (compared to control) will increase sharing intentions

✅ H5b: Other-focused intervention (compared to control) will increase sharing intentions

We replicate our previous work in the comment group: the self- and other-focused conditions increase sharing intentions compared to the control, and these effects are stronger for narrowcast compared to broadcasting sharing intentions

✅ These effects were smaller for both the timed and untimed groups

mod_h5 = lmer(msg_share ~ 1 + article_cond*condition + (1 + article_cond | SID),
              data = data,
              control = lmerControl(optimizer = "bobyqa"))

plot

predicted_h5 = ggeffects::ggpredict(mod_h5, c("article_cond", "condition")) %>%
              data.frame() %>%
  mutate(model = "sharing")

predicted_h5 %>%
  mutate(x = factor(x, levels = c("self", "control", "other"))) %>%
  ggplot(aes(x = group, y = predicted, color = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1) +
  coord_flip() +
  scale_color_manual(name = "", values = palette_condition) +
  labs(x = "", y = "\nmean predicted sharing intention rating") +
  plot_aes +
  theme(legend.position = "top")

model table

table_h5 = table_model(mod_h5)

table_h5 %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 32.83 [28.09, 37.56] 451.92 13.62 < .001
other 14.64 [11.26, 18.02] 460.57 8.51 < .001
self 9.46 [6.16, 12.76] 460.76 5.64 < .001
group (timed) 0.70 [-5.61, 7.02] 447.26 0.22 .827
group (untimed) -0.20 [-6.43, 6.04] 447.97 -0.06 .951
other x group (timed) -13.68 [-18.16, -9.20] 449.95 -6.00 < .001
self x group (timed) -8.41 [-12.79, -4.04] 451.59 -3.78 < .001
other x group (untimed) -13.32 [-17.74, -8.89] 450.29 -5.91 < .001
self x group (untimed) -9.26 [-13.59, -4.94] 452.15 -4.21 < .001

simple effects by group

simple_effects(mod_h5, sharing = FALSE)
contrast group estimate p.value
other - control comment 14.64 < .001
self - control comment 9.46 < .001
other - control timed 0.96 .522
self - control timed 1.05 .475
other - control untimed 1.32 .364
self - control untimed 0.20 .890

summary

summary(mod_h5)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_share ~ 1 + article_cond * condition + (1 + article_cond |  
##     SID)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 49609.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4244 -0.5058 -0.0990  0.5157  3.7757 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr       
##  SID      (Intercept)       584.02   24.166              
##           article_condother 100.21   10.011    0.01      
##           article_condself   84.29    9.181   -0.08  0.28
##  Residual                   503.50   22.439              
## Number of obs: 5312, groups:  SID, 448
## 
## Fixed effects:
##                                    Estimate Std. Error       df t value
## (Intercept)                         32.8256     2.4099 451.9229  13.621
## article_condother                   14.6377     1.7191 460.5726   8.515
## article_condself                     9.4588     1.6781 460.7627   5.637
## conditiontimed                       0.7031     3.2139 447.2603   0.219
## conditionuntimed                    -0.1958     3.1732 447.9713  -0.062
## article_condother:conditiontimed   -13.6791     2.2797 449.9517  -6.000
## article_condself:conditiontimed     -8.4126     2.2263 451.5892  -3.779
## article_condother:conditionuntimed -13.3151     2.2529 450.2884  -5.910
## article_condself:conditionuntimed   -9.2624     2.1997 452.1467  -4.211
##                                                Pr(>|t|)    
## (Intercept)                        < 0.0000000000000002 ***
## article_condother                  0.000000000000000238 ***
## article_condself                   0.000000030263850019 ***
## conditiontimed                                 0.826939    
## conditionuntimed                               0.950814    
## article_condother:conditiontimed   0.000000004052301585 ***
## article_condself:conditiontimed                0.000179 ***
## article_condother:conditionuntimed 0.000000006750561056 ***
## article_condself:conditionuntimed  0.000030724212555648 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                     (Intr) artcl_cndt artcl_cnds cndtnt cndtnn
## artcl_cndth         -0.258                                    
## artcl_cndsl         -0.305  0.451                             
## conditintmd         -0.750  0.193      0.228                  
## conditnntmd         -0.759  0.196      0.231      0.569       
## artcl_cndthr:cndtnt  0.194 -0.754     -0.340     -0.252 -0.148
## artcl_cndslf:cndtnt  0.230 -0.340     -0.754     -0.300 -0.174
## artcl_cndthr:cndtnn  0.197 -0.763     -0.344     -0.148 -0.253
## artcl_cndslf:cndtnn  0.232 -0.344     -0.763     -0.174 -0.301
##                     artcl_cndthr:cndtnt artcl_cndslf:cndtnt artcl_cndthr:cndtnn
## artcl_cndth                                                                    
## artcl_cndsl                                                                    
## conditintmd                                                                    
## conditnntmd                                                                    
## artcl_cndthr:cndtnt                                                            
## artcl_cndslf:cndtnt  0.447                                                     
## artcl_cndthr:cndtnn  0.575               0.259                                 
## artcl_cndslf:cndtnn  0.259               0.575               0.447

combined plot

versus control

predicted_h2a_con = ggeffects::hypothesis_test(mod_h2a, c("article_cond", "condition")) %>%
  data.frame() %>%
  arrange(condition) %>%
  filter(grepl("^control", article_cond) & !article_cond == "control-control") %>%
  filter(condition %in% c("comment-comment", "untimed-untimed", "timed-timed")) %>%
  mutate(model = "self-relevance")

predicted_h2b_con = ggeffects::hypothesis_test(mod_h2b, c("article_cond", "condition")) %>%
  data.frame() %>%
  arrange(condition) %>%
  filter(grepl("^control", article_cond) & !article_cond == "control-control") %>%
  filter(condition %in% c("comment-comment", "untimed-untimed", "timed-timed")) %>%
  mutate(model = "social relevance")

predicted_h5_con = ggeffects::hypothesis_test(mod_h5, c("article_cond", "condition")) %>%
  data.frame() %>%
  arrange(condition) %>%
  filter(grepl("^control", article_cond) & !article_cond == "control-control") %>%
  filter(condition %in% c("comment-comment", "untimed-untimed", "timed-timed")) %>%
  mutate(model = "sharing")

bind_rows(predicted_h2a_con, predicted_h2b_con, predicted_h5_con) %>%
  mutate(condition = gsub("-.*", "", condition),
         condition = recode(condition, "untimed" = "reflect:\nuntimed",
                            "timed" = "reflect:\ntimed"),
         article_cond = gsub("control-", "", article_cond),
         article_cond = sprintf("%s > control", article_cond),
         Contrast = Contrast * -1,
         conf.low = conf.low * -1,
         conf.high = conf.high * -1,
         article_cond = factor(article_cond, levels = c("self > control", "other > control"))) %>%
  mutate(model = factor(model, levels = c("self-relevance", "social relevance", "sharing"))) %>%
  ggplot(aes(x = condition, y = Contrast, color = article_cond)) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = .75) +
  facet_grid(~ model) +
  coord_flip() +
  scale_color_manual(name = "", values = c(palette_condition[[1]], palette_condition[[3]])) +
  labs(x = "", y = "\npredicted rating") +
  plot_aes +
  theme(legend.position = "top")

all conditions

bind_rows(predicted_h2a, predicted_h2b, predicted_h5) %>%
  mutate(model = factor(model, levels = c("self-relevance", "social relevance", "sharing")),
         x = factor(x, levels = c("self", "control", "other")),
         group = ifelse(group == "timed", "reflect:\ntimed",
                 ifelse(group == "untimed", "reflect:\nuntimed", "comment")),
         group = factor(group, levels = c("reflect:\ntimed", "reflect:\nuntimed", "comment"))) %>%
  ggplot(aes(x = group, y = predicted, color = x)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), position = position_dodge(.5), size = 1) +
  facet_grid(~ model) +
  coord_flip() +
  scale_color_manual(name = "", values = palette_condition) +
  labs(x = "", y = "\npredicted rating") +
  plot_aes +
  theme(legend.position = "top")

combined table

table_h2a %>% mutate(DV = "H2a: Self-relevance") %>%
  bind_rows(table_h2b %>% mutate(DV = "H2b: Social relevance")) %>%
  bind_rows(table_h5 %>% mutate(DV = "H5: Sharing intention")) %>%
  select(DV, everything()) %>%
  kable() %>%
  kable_styling()
DV term b [95% CI] df t p
H2a: Self-relevance intercept 49.91 [45.62, 54.20] 456.24 22.85 < .001
H2a: Self-relevance other 4.03 [0.43, 7.62] 462.82 2.20 .028
H2a: Self-relevance self 14.44 [10.47, 18.40] 457.92 7.15 < .001
H2a: Self-relevance group (timed) -3.52 [-9.23, 2.19] 449.41 -1.21 .227
H2a: Self-relevance group (untimed) -4.66 [-10.31, 0.98] 450.11 -1.62 .105
H2a: Self-relevance other x group (timed) -5.68 [-10.44, -0.92] 451.05 -2.35 .019
H2a: Self-relevance self x group (timed) -9.83 [-15.10, -4.56] 448.78 -3.67 < .001
H2a: Self-relevance other x group (untimed) -3.26 [-7.96, 1.45] 451.42 -1.36 .174
H2a: Self-relevance self x group (untimed) -11.66 [-16.86, -6.46] 449.43 -4.40 < .001
H2b: Social relevance intercept 53.76 [49.40, 58.12] 455.12 24.22 < .001
H2b: Social relevance other 14.45 [10.91, 17.99] 459.95 8.03 < .001
H2b: Social relevance self 9.87 [6.28, 13.47] 459.93 5.39 < .001
H2b: Social relevance group (timed) -3.11 [-8.92, 2.70] 449.05 -1.05 .293
H2b: Social relevance group (untimed) -4.36 [-10.10, 1.38] 449.79 -1.49 .136
H2b: Social relevance other x group (timed) -13.98 [-18.67, -9.29] 448.90 -5.86 < .001
H2b: Social relevance self x group (timed) -6.86 [-11.63, -2.09] 450.70 -2.83 .005
H2b: Social relevance other x group (untimed) -10.26 [-14.89, -5.62] 449.24 -4.35 < .001
H2b: Social relevance self x group (untimed) -7.39 [-12.11, -2.68] 451.27 -3.08 .002
H5: Sharing intention intercept 32.83 [28.09, 37.56] 451.92 13.62 < .001
H5: Sharing intention other 14.64 [11.26, 18.02] 460.57 8.51 < .001
H5: Sharing intention self 9.46 [6.16, 12.76] 460.76 5.64 < .001
H5: Sharing intention group (timed) 0.70 [-5.61, 7.02] 447.26 0.22 .827
H5: Sharing intention group (untimed) -0.20 [-6.43, 6.04] 447.97 -0.06 .951
H5: Sharing intention other x group (timed) -13.68 [-18.16, -9.20] 449.95 -6.00 < .001
H5: Sharing intention self x group (timed) -8.41 [-12.79, -4.04] 451.59 -3.78 < .001
H5: Sharing intention other x group (untimed) -13.32 [-17.74, -8.89] 450.29 -5.91 < .001
H5: Sharing intention self x group (untimed) -9.26 [-13.59, -4.94] 452.15 -4.21 < .001

word count analyses

Test whether word count is higher in the experimental conditions, and whether it is positively associated with self and social relevance, and sharing intention ratings.

descriptives

words_ratings = data %>%
  filter(condition == "comment") %>%
  ungroup() %>%
  mutate(n_c = n_words - mean(n_words, na.rm = TRUE))

data %>%
  group_by(article_cond) %>%
  summarize(mean = mean(n_words, na.rm = TRUE),
            sd = sd(n_words, na.rm = TRUE),
            min = min(n_words, na.rm = TRUE),
            max = max(n_words, na.rm = TRUE)) %>%
  kable(digits = 2) %>%
  kableExtra::kable_styling()
article_cond mean sd min max
control 14.23 7.21 6 72
other 17.62 9.22 6 69
self 18.47 10.33 6 72

condition effects

Is word count higher in the experimental conditions compared to the control condition?

✅ The word count is higher in the experimental conditions compared to the control condition

mod_words = lmer(n_words ~ 1 + article_cond + (1 + article_cond | SID),
              data = data,
              control = lmerControl(optimizer = "bobyqa"))

plot

predicted_words = ggeffects::ggpredict(mod_words, c("article_cond")) %>%
              data.frame() %>%
  mutate(x = factor(x, levels = c("self", "control", "other")))

predicted_sub_words = ggeffects::ggpredict(mod_words, terms = c("article_cond", "SID"), type = "random") %>%
  data.frame()%>%
  mutate(x = factor(x, levels = c("self", "control", "other")))

predicted_words %>%
  ggplot(aes(x = x, y = predicted)) +
  stat_summary(data = predicted_sub_words, aes(group = group), fun = "mean", geom = "line",
               size = .08, color = "grey50") +
  stat_summary(aes(group = group), fun = "mean", geom = "line", size = 1) +
  geom_pointrange(aes(color = x, ymin = conf.low, ymax = conf.high), size = .5) +
  scale_color_manual(name = "", values = palette_condition, guide = "none") +
  scale_alpha_manual(name = "", values = c(1, .5)) +
  labs(x = "", y = "predicted word count\n") +
  plot_aes

model table

table_words = table_model(mod_words, word_count = TRUE)

table_words %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 14.02 [13.01, 15.03] 122.32 27.48 < .001
other 3.38 [2.29, 4.47] 123.52 6.15 < .001
self 4.23 [3.01, 5.44] 124.00 6.87 < .001

summary

summary(mod_words)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: n_words ~ 1 + article_cond + (1 + article_cond | SID)
##    Data: data
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 9627.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8325 -0.5108 -0.1123  0.3834  6.5006 
## 
## Random effects:
##  Groups   Name              Variance Std.Dev. Corr     
##  SID      (Intercept)       22.90    4.786             
##           article_condother 19.02    4.361    0.19     
##           article_condself  28.61    5.348    0.34 0.72
##  Residual                   34.85    5.903             
## Number of obs: 1434, groups:  SID, 124
## 
## Fixed effects:
##                   Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)        14.0201     0.5103 122.3236  27.477 < 0.0000000000000002 ***
## article_condother   3.3815     0.5495 123.5169   6.154        0.00000000970 ***
## article_condself    4.2266     0.6155 124.0030   6.867        0.00000000028 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) artcl_cndt
## artcl_cndth -0.157           
## artcl_cndsl -0.018  0.621

relevance

Is word count positively associated with self and social relevance ratings?

self-relevance

✅ Word count is positively associated with self-relevance ratings

mod_words_h1 = lmer(msg_rel_self ~ 1 + n_c + (1 + n_c | SID),
              data = filter(words_ratings, sharing_type == 0),
              control = lmerControl(optimizer = "bobyqa"))

plot

values = seq(-15, 60, 10)
predicted_self = ggeffects::ggpredict(mod_words_h1, terms = "n_c [values]") %>%
  data.frame()

predicted_sub_self = ggeffects::ggpredict(mod_words_h1, terms = c("n_c [values]", "SID"), type = "random") %>%
  data.frame()

predicted_self %>%
  ggplot(aes(x, predicted)) +
  stat_smooth(data = predicted_sub_self, aes(group = group), geom ='line', method = "lm", alpha = .05, size = 1, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5) +
  geom_line(size = 1) +
  labs(x = "\nword count (grand mean-centered)", y = "predicted self-relevance rating\n") +
  plot_aes

model table

table_words_h1 = table_model(mod_words_h1, word_count = TRUE)

table_words_h1 %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 56.21 [52.13, 60.29] 123.02 27.27 < .001
word count 0.49 [0.23, 0.74] 56.53 3.81 < .001

summary

summary(mod_words_h1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_self ~ 1 + n_c + (1 + n_c | SID)
##    Data: filter(words_ratings, sharing_type == 0)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 13779.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9261 -0.6255  0.1300  0.6758  2.7583 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  SID      (Intercept) 452.1796 21.2645      
##           n_c           0.3112  0.5578  0.18
##  Residual             722.5498 26.8803      
## Number of obs: 1433, groups:  SID, 124
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  56.2083     2.0611 123.0217  27.271 < 0.0000000000000002 ***
## n_c           0.4864     0.1277  56.5310   3.809             0.000346 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## n_c 0.122

social relevance

✅ Word count is positively associated with social relevance ratings

mod_words_h2 = lmer(msg_rel_social ~ 1 + n_c + (1 + n_c | SID),
              data = filter(words_ratings, sharing_type == 0),
              control = lmerControl(optimizer = "bobyqa"))

plot

values = seq(-15, 60, 10)
predicted_social = ggeffects::ggpredict(mod_words_h2, terms = "n_c [values]") %>%
  data.frame()

predicted_sub_social = ggeffects::ggpredict(mod_words_h2, terms = c("n_c [values]", "SID"), type = "random") %>%
  data.frame()

predicted_social %>%
  ggplot(aes(x, predicted)) +
  stat_smooth(data = predicted_sub_social, aes(group = group), geom ='line', method = "lm", alpha = .05, size = 1, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5) +
  geom_line(size = 1) +
  labs(x = "\nword count (grand mean-centered)", y = "predicted social relevance rating\n") +
  plot_aes

model table

table_words_h2 = table_model(mod_words_h2, word_count = TRUE)

table_words_h2 %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 62.13 [57.96, 66.31] 123.19 29.47 < .001
word count 0.46 [0.23, 0.69] 76.31 3.95 < .001

summary

summary(mod_words_h2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_rel_social ~ 1 + n_c + (1 + n_c | SID)
##    Data: filter(words_ratings, sharing_type == 0)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 13459
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0993 -0.4577  0.1436  0.6199  2.8391 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  SID      (Intercept) 491.7477 22.1754      
##           n_c           0.2901  0.5386  0.20
##  Residual             561.9133 23.7047      
## Number of obs: 1433, groups:  SID, 124
## 
## Fixed effects:
##             Estimate Std. Error       df t value             Pr(>|t|)    
## (Intercept)  62.1337     2.1085 123.1917  29.468 < 0.0000000000000002 ***
## n_c           0.4615     0.1169  76.3134   3.949             0.000173 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## n_c 0.136

sharing

Is word count positively associated with sharing intention ratings?

Here we focus on narrowcasting only since that is the type of sharing we used in fMRI study 1.

✅ Word count is positively associated with narrowcast sharing intentions

mod_words_h3 = lmer(msg_share ~ 1 + n_c + (1 + n_c | SID),
              data = filter(words_ratings, sharing_type == 0),
              control = lmerControl(optimizer = "bobyqa"))

plot

values = seq(-15, 60, 10)
predicted_sharing = ggeffects::ggpredict(mod_words_h3, terms = "n_c [values]") %>%
  data.frame()

predicted_sub_sharing = ggeffects::ggpredict(mod_words_h3, terms = c("n_c [values]", "SID"), type = "random") %>%
  data.frame()

predicted_sharing %>%
  ggplot(aes(x, predicted)) +
  stat_smooth(data = predicted_sub_sharing, aes(group = group), geom ='line', method = "lm", alpha = .05, size = 1, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .5) +
  geom_line(size = 1) +
  labs(x = "\nword count (grand mean-centered)", y = "predicted sharing intention rating\n") +
  plot_aes

model table

table_words_h3 = table_model(mod_words_h3, word_count = TRUE)

table_words_h3 %>%
    kable()  %>%
    kableExtra::kable_styling()
term b [95% CI] df t p
intercept 41.25 [36.52, 45.97] 122.59 17.28 < .001
word count 0.33 [0.07, 0.58] 49.54 2.54 .014

summary

summary(mod_words_h3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: msg_share ~ 1 + n_c + (1 + n_c | SID)
##    Data: filter(words_ratings, sharing_type == 0)
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 13559.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.00445 -0.56241 -0.09983  0.54949  2.85174 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr
##  SID      (Intercept) 640.3267 25.3047      
##           n_c           0.4505  0.6712  0.15
##  Residual             588.6676 24.2625      
## Number of obs: 1433, groups:  SID, 124
## 
## Fixed effects:
##             Estimate Std. Error       df t value            Pr(>|t|)    
## (Intercept)  41.2461     2.3873 122.5948  17.277 <0.0000000000000002 ***
## n_c           0.3252     0.1279  49.5373   2.542              0.0142 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##     (Intr)
## n_c 0.128

combined plot

predicted_sub = predicted_sub_self %>%
  mutate(facet = "self-relevance") %>%
  bind_rows(., predicted_sub_social %>%  mutate(facet = "social relevance")) %>%
  bind_rows(., predicted_sub_sharing %>%  mutate(facet = "sharing")) %>%
  mutate(facet = factor(facet, levels = c("self-relevance", "social relevance", "sharing")))
  
predicted_self %>%
  mutate(facet = "self-relevance") %>%
  bind_rows(., predicted_social %>%  mutate(facet = "social relevance")) %>%
  bind_rows(., predicted_sharing %>%  mutate(facet = "sharing")) %>%
  mutate(facet = factor(facet, levels = c("self-relevance", "social relevance", "sharing"))) %>%
  ggplot(aes(x, predicted, color = facet, fill = facet)) +
  stat_smooth(data = predicted_sub, aes(group = interaction(group, facet)), geom ='line', method = "lm", alpha = .1, size = .75, se = FALSE) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25, color = NA) +
  geom_line(size = 2) +
  scale_color_manual(values = palette_dv, name = "") + 
  scale_fill_manual(values = palette_dv, name = "") + 
  labs(x = "\nword count (grand mean-centered)", y = "predicted rating\n") +
  plot_aes +
  theme(legend.position = c(.2, .95))

combined table

table_words %>% mutate(DV = "Word count") %>%
  bind_rows(table_words_h1 %>% mutate(DV = "Self-relevance")) %>%
  bind_rows(table_words_h2 %>% mutate(DV = "Social relevance")) %>%
  bind_rows(table_words_h3 %>% mutate(DV = "Sharing intention")) %>%
  select(DV, everything()) %>%
  kable() %>%
  kable_styling()
DV term b [95% CI] df t p
Word count intercept 14.02 [13.01, 15.03] 122.32 27.48 < .001
Word count other 3.38 [2.29, 4.47] 123.52 6.15 < .001
Word count self 4.23 [3.01, 5.44] 124.00 6.87 < .001
Self-relevance intercept 56.21 [52.13, 60.29] 123.02 27.27 < .001
Self-relevance word count 0.49 [0.23, 0.74] 56.53 3.81 < .001
Social relevance intercept 62.13 [57.96, 66.31] 123.19 29.47 < .001
Social relevance word count 0.46 [0.23, 0.69] 76.31 3.95 < .001
Sharing intention intercept 41.25 [36.52, 45.97] 122.59 17.28 < .001
Sharing intention word count 0.33 [0.07, 0.58] 49.54 2.54 .014

cite packages

report::cite_packages()
##   - Angelo Canty and Brian Ripley (2021). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-28.
##   - Douglas Bates, Martin Maechler and Mikael Jagan (2023). Matrix: Sparse and Dense Matrix Classes and Methods. R package version 1.5-4. https://CRAN.R-project.org/package=Matrix
##   - Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
##   - Evan Kleiman (2021). EMAtools: Data Management Tools for Real-Time Monitoring/Ecological Momentary Assessment Data. R package version 0.1.4. https://CRAN.R-project.org/package=EMAtools
##   - H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
##   - Hadley Wickham (2021). forcats: Tools for Working with Categorical Variables (Factors). R package version 0.5.1. https://CRAN.R-project.org/package=forcats
##   - Hadley Wickham (2022). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.0. https://CRAN.R-project.org/package=stringr
##   - Hadley Wickham and Maximilian Girlich (2022). tidyr: Tidy Messy Data. R package version 1.2.0. https://CRAN.R-project.org/package=tidyr
##   - Hadley Wickham, Jennifer Bryan and Malcolm Barrett (2021). usethis: Automate Package and Project Setup. R package version 2.1.5. https://CRAN.R-project.org/package=usethis
##   - Hadley Wickham, Jim Hester and Jennifer Bryan (2022). readr: Read Rectangular Text Data. R package version 2.1.2. https://CRAN.R-project.org/package=readr
##   - Hadley Wickham, Jim Hester, Winston Chang and Jennifer Bryan (2021). devtools: Tools to Make Developing R Packages Easier. R package version 2.4.3. https://CRAN.R-project.org/package=devtools
##   - Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2022). dplyr: A Grammar of Data Manipulation. R package version 1.0.9. https://CRAN.R-project.org/package=dplyr
##   - Hao Zhu (2021). kableExtra: Construct Complex Table with 'kable' and Pipe Syntax. R package version 1.3.4. https://CRAN.R-project.org/package=kableExtra
##   - Jim Hester, Hadley Wickham and Gábor Csárdi (2021). fs: Cross-Platform File System Operations Based on 'libuv'. R package version 1.5.2. https://CRAN.R-project.org/package=fs
##   - Kirill Müller and Hadley Wickham (2022). tibble: Simple Data Frames. R package version 3.1.8. https://CRAN.R-project.org/package=tibble
##   - Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package:Tests in Linear Mixed Effects Models." _Journal of StatisticalSoftware_, *82*(13), 1-26. doi: 10.18637/jss.v082.i13 (URL:https://doi.org/10.18637/jss.v082.i13).
##   - Lionel Henry and Hadley Wickham (2020). purrr: Functional Programming Tools. R package version 0.3.4. https://CRAN.R-project.org/package=purrr
##   - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects fromRegression Models." _Journal of Open Source Software_, *3*(26), 772.doi: 10.21105/joss.00772 (URL: https://doi.org/10.21105/joss.00772).
##   - Lüdecke D, Ben-Shachar M, Patil I, Makowski D (2020). "Extracting,Computing and Exploring the Parameters of Statistical Models using R."_Journal of Open Source Software_, *5*(53), 2445. doi:10.21105/joss.02445 (URL: https://doi.org/10.21105/joss.02445).
##   - R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
##   - Rinker, T. W. & Kurkiewicz, D. (2017). pacman: Package Management for R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman
##   - Wickham et al., (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686
##   - Yihui Xie (2021). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.37.